% Mutual information between photon number difference and phase difference


clear

nmax=50;
dfi=0.05;
ns=-nmax:nmax;
fis=-pi:dfi:pi;
nsiz=size(ns,2);
fisiz=size(fis,2);

sigman=3;

[N,FI]=meshgrid(ns,fis);

probs0=exp(-(N.^2/2/sigman^2));
probs0=probs0/sum(sum(probs0))/dfi;
probsN0=sum(probs0)*dfi;
probsfi0=sum(probs0')';
entrN0=-sum(probsN0.*log(probsN0))
entrfi0=-sum(probsfi0.*log(probsfi0))*dfi
celkentr0=entrN0+entrfi0
entrtot0=-sum(sum(probs0.*log(probs0)))*dfi



ampln=10.6;
probs1=exp(-((N-ampln*cos(FI)).^2/2/sigman^2));
probs1=probs1/sum(sum(probs1))/dfi;


probsN1=sum(probs1)*dfi;
probsfi1=sum(probs1')';
productprobs1=probsfi1*probsN1;

%size(productprobs1)
%sum(sum(productprobs1))*dfi

logprobs1=log(probs1)-log(productprobs1);
inform1=sum(sum(probs1.*logprobs1))*dfi
entrN1=-sum(probsN1.*log(probsN1))
entrfi1=-sum(probsfi1.*log(probsfi1))*dfi
celkentr1=entrN1+entrfi1-inform1
entrtot1=-sum(sum(probs1.*log(probs1)))*dfi


% Posun Kerrem
chi=0.13;
probs2=zeros(size(probs1));
for n=1:nsiz
 shift=chi*ns(n);
 discrshift=round(shift/dfi);
 probs2(:,n)=circshift(probs1(:,n),discrshift);
end

probsN2=sum(probs2)*dfi;
probsfi2=sum(probs2')';
productprobs2=probsfi2*probsN2;
logprobs2=log(probs2)-log(productprobs2);
inform2=sum(sum(probs2.*logprobs2))*dfi
entrN2=-sum(probsN2.*log(probsN2))
entrfi2=-sum(probsfi2.*log(probsfi2))*dfi
celkentr2=entrN2+entrfi2-inform2
entrtot2=-sum(sum(probs2.*log(probs2)))*dfi

ninmean=500;
nouts=0:2.5*ninmean;
noutsiz=size(nouts,2);
probsnout=zeros(1,noutsiz);
probsnoutp=zeros(1,noutsiz);
probsnoutmin=zeros(1,noutsiz);

for n=1:fisiz
 pfi=fis(n);
 npmean=ninmean*(1+sin(pfi))+1;
 sigp=sqrt(npmean);
 pomprob=exp(-(nouts-npmean).^2/2/npmean);
 pomprob=pomprob/sum(pomprob);
 probsnout=probsnout+probsfi2(n)*pomprob;
 probsnoutp=probsnoutp+probsfi1(n)*pomprob;
 npmean=ninmean*(1-sin(pfi))+1;
 sigp=sqrt(npmean);
 pomprob=exp(-(nouts-npmean).^2/2/npmean);
 pomprob=pomprob/sum(pomprob);
 probsnoutmin=probsnoutmin+probsfi2(n)*pomprob;
end
probsnout=probsnout/sum(probsnout);
meannout=probsnout*nouts';
probsnoutp=probsnoutp/sum(probsnoutp);
meannoutp=probsnoutp*nouts';
probsnoutmin=probsnoutmin/sum(probsnoutmin);
meannoutmin=probsnoutmin*nouts';

velik=20;

figure(1);
%contour(fis,ns,probs1')
mesh(fis,ns,probs1')
set(gca,'FontSize',velik)
%xlabel('\Delta \phi','FontSize',velik)
%ylabel('\Delta n','FontSize',velik)
axis([-pi,pi,-nmax,nmax,0,.025])
view([-20,50])

figure(2)
plot(ns,[probsN0;probsN1;probsN2])
xlabel('\Delta n','FontSize',velik)
ylabel('probability','FontSize',velik)
set(gca,'FontSize',velik)

figure(3);
mesh(fis,ns,probs2')
%surf(fis,ns,probs2')
set(gca,'FontSize',velik)
axis([-pi,pi,-nmax,nmax,0,.025])
view([-20,50])

figure(4)
plot(fis,probsfi1','r')
hold on
plot(fis,probsfi2','b')
hold off
axis([-pi,pi,0,.6])
%xlabel('\Delta \phi','FontSize',velik)
%ylabel('probability density','FontSize',velik)
set(gca,'FontSize',velik)

figure(5)
plot(nouts,probsnout,'b')
hold on
plot([meannout,meannout],[0,.005],'b--')
plot(nouts,probsnoutp,'r')
plot([meannoutp,meannoutp],[0,.005],'r--')
plot(nouts,probsnoutmin,'g')
plot([meannoutmin,meannoutmin],[0,.005],'g--')
hold off
axis([0,2.2*ninmean,0,.01])
set(gca,'FontSize',velik)


